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We present a new approach to find accurate solutions to the Poisson equation, as ob¬ 
tained from the steady-state limit of a diffusion equation with strong source terms. For 
this purpose, we start from Boltzmann’s kinetic theory and investigate the influence of 
higher order terms on the resulting macroscopic equations. By performing an appropriate 
expansion of the equilibrium distribution, we provide a method to remove the unneces¬ 
sary terms up to a desired order and show that it is possible to find, with high level of 
accuracy, the steady-state solution of the diffusion equation for sizeable Knudsen num¬ 
bers. In order to test our kinetic approach, we discretise the Boltzmann equation and 
solve the Poisson equation, spending up to six order of magnitude less computational 
time for a given precision than standard lattice Boltzmann methods. 

Keywords: High Knudsen number, higher-order moments, diffusion equation, Poisson 
equation, lattice Boltzmann 


1. Introduction 


The diffusion equation is widely used to describe transport phenomena in which 
heat, mass, and other physical quantities are transferred in space and time due to 
underlying molecular collision processes &. At steady-state and when the source 
term does not depend explicitly on the concentration of the associated field, such 
diffusion processes settle down in the form of a non-local relation between the spatial 
distribution of the source and the resulting concentration of the associated field, a 
relation which is governed by the Poisson equation The notion of the Poisson 
equation as the steady-state solution of the diffusion equation extends to many other 
phenomena, e.g. electrostatics systems plasma physics SI, chemistry and biology 
molecular biophysics^, density functional theory and solid-state physicsEMSII, 
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to name but a few. In fact, solving the Poisson equation accurately and efficiently 
continues to be a subject of intense research to this day. 

From a microscopic point of view, the diffusion equation emerges from the un¬ 
derlying Boltzmann kinetic equation in the limit of very small Knudsen numbers, 
i.e. molecular mean free path much shorter than the typical macroscopic scale. Solv¬ 
ing the Boltzmann kinetic equation to study diffusion processes makes apparently 
little sense, since the Boltzmann distribution function lives in a double-dimensional 
phase-space, i.e. position and velocity. However, in the last two decades, minimal 
(lattice) versions of the Boltzmann equation have been developed, in which veloc¬ 
ity space is reduced to a small set of discrete velocities, so that the computational 
cost is cut down dramatically, making the solution of the lattice kinetic equation 
often competitive towards standard grid-discretisation of the corresponding partial 
differential equation. 

The Lattice Boltzmann Method (LBM) has attracted much interest for solving 
the Navier-Stokes equations and h as been extended to solve other type of 

phenomena, e.g. diffusion equation USDS Maxwell equation J22 quantum systems 
■22, relativistic fluid dynamics 2? I 18 I, and the Poisson equationE ^9 | 20 | 21 | 22 |, name 
but a few. To the best of our knowledge, to date, the solution of the Poisson equation 
based on kinetic theory ha s been confined to smooth (non-stiff) source terms and 
small Knudsen numbers ^3. if the target is the steady-state solution of the diffusion 
equation, as it is the case for the Poisson equation, the kinetic approach looses 
much of its appeal, since its inherently real-time dynamics must proceed through 
smaller time-steps than those affordable by fictitious-time iteration techniques. In 
this paper we shall show that such gap can be closed by resorting to higher-order 
kinetic formulations, at least for the cases where the source term does not depend 
on the concentration field avoiding artificial solutions as pointed out in Ref. 

Higher-order formulations of the LB equation have been developed before, 
mostly in the context of thermal ^and relativistic fluid dynamics ^51, an( f recently, 
density functional theory 23 There, the main idea is to design the equilibrium dis¬ 
tribution in such a way to include not only the conserved moments, such as density 
and current, but also the non-conserved ones, namely the momentum flux tensor 
and the heat flux, and eventually higher order kinetic moments with no direct hy¬ 
drodynamic significance (sometimes called “ghosts”). Clearly, in order to match 
this longer ladder of kinetic moments, a correspondingly larger set of discrete ve¬ 
locities is required. It should be emphasized that the discrete kinetic equation is a 
superset of the diffusion equation, and consequently, to the purpose of solving the 
diffusion equation alone, higher-order contributions must be regarded as discretisa¬ 
tion artefacts which need to be minimised, and possibly canceled altogether. This is 
precisely the target of this paper. Here, we will use higher-order lattices and associ¬ 
ated equilibria, to remove higher-order contributions to the macroscopic equations 
reproduced by the Boltzmann equation. For this purpose, we develop the respective 
LBM and show that, by expanding the equilibrium distribution up to seventh order 
in Hermite polynomials, one can solve the Poisson equation six orders of magnitude 
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faster than with standard LB models, at a given level accuracy. 

The paper is divided as follows: In Sec. [2j we introduce the diffusion equation 
with stiff source term and the associated Boltzmann kinetic equation. In Sec. [3l we 
describe the influence of higher order moments on the steady-state solution of the 
diffusion equation. In Sec. [2 we develop a LBM for the Poisson equation, validate 
and compare the model with standard cases. Finally, we discuss the results and 
outlook of our work. 


2. From Boltzmann kinetic theory to diffusion equation 

Let us begin by writing the diffusion equation for a scalar $ in the standard form: 

^=DVH + S , ( 1 ) 


where $ is the concentration of particles (temperature in the case of the heat trans¬ 
fer equation), D is the diffusivity, and S is the source term, which, in our case, 
depending neither on $ nor on time. The steady-state version of Eq. © delivers 
the Poisson equation with source —S/D. It is well-known that this equation can be 
obtained from an underlying Boltzmann kinetic equation, via a Chapman-Enskog 
expansion, in the limit of small Knudsen number s, K n = X/L, being A the mean-free 
path and L the characteristic size of the system ^ A 

The Boltzmann equation in the single relaxation time approximation (BGK) ^H1 
reads as follows: 

% + v.Vf = - 1 -(f-r) + S , (2) 

where / = f(x,v,t) is the probability distribution function, v are the microscopic 
velocity vectors, r is the relaxation time, and / eq is the equilibrium distribution. 

In the non-relativistic context, the local equilibrium is given by a Maxwell- 
Boltzmann distribution: 


/ 


eq _ 


e -{v-u) 2 /2 c 2 

(2^ C 2)3/2 


( 3 ) 


where u is the local flow speed. Such local equilibrium encodes the basic symmetries 
of Newtonian mechanics, particularly Galilean invariance, which is built-in via the 
dependence on the molecular speed relative to the fluid v — u, rather than the 
absolute one, v. In order to secure Galilean invariance for any fluid speed, kinetic 
moments at all orders should be matched on the lattice. This would require an 
infinite connectivity, which is clearly unviable on any realistic lattice. As a result, 
LB l oca l equilibria are typically based on finite-order Hermite expansions, of the 
form [21 


r(x, v, t ) = w(v) £ a,™ (x, t)H <"> (v) , (4) 

n —0 
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where Hr, 1 ' 1 are the tensorial Hermite polynomials of order n, and w(v) is the weight 
function defined as, 


w(v) 


(2nc 2 s ) 3 / 2 


0 -v 2 /2c 2 


(5) 


c s being the sound speed. Note that w(v) is the uniform global equilibrium, corre¬ 
sponding the no-flow limit of the local equilibrium / eq . The coefficients (kinetic 
moments) are calculated by projecting the equilibrium distribution onto the respec¬ 
tive Hermite polynomial, 


a (n) = I f e qH (n) d 3 y ( 6 ) 

The lowest-order kinetic moments have a direct macroscopic interpretation. For 
instance, by integrating on the velocity space (Ho = 1), 


$ = 



/ 


/ eq d 3 v 


(7) 


where the last equality comes from the requirement of mass conservation on the 
collsion operator, namely: 

\ /(/ - D d 3 v = 0 . (8) 

The diffusivity of the system is dictated by the relaxation time as 

D = c 2 s t. (9) 


It can be shown that the Boltzmann equation recovers the diffusion equation 
in the limit of small Knudsen numbers, which is related with the mean-free path 
and therefore with the relaxation time. This is the strongly-interacting fluid regime, 
whereby the particle dynamics is collision-dominated, so that collective behaviour 
sets in on a short timescale r. The opposite limit of large Knudsen numbers, hence 
large values of r, corresponds to the free-particle motion (ballistic regime), whereby 
memory of the initial conditions is kept for a very long time. It is therefore clear that 
in the ballistic regime higher order moments of the equilibrium distribution do not 
relax and consequently the Boltzmann equation does not converge to any standard 
diffusion equation. On the other hand, by increasing the relaxation time, hence the 
effective diffusivity, one would intuitively expect a quickest path to steady-state. In 
this respect, it would be highly desirable to derive a diffusion equation in kinetic 
form, which can attain steady-state solution at relatively large Knudsen numbers 
by exactly cancelling as many high-order terms as possible. 

As discussed above, this can be achieved by zeroing as many higher order con¬ 
tributions as possible, while retaining the largest possible diffusion coefficient. 
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3. Influence of Higher-Order Terms 

Let us consider the case when all time derivatives are zero, d/dt = 0. The Boltzmann 
equation, Eq. m, can be written as, 

f = n + TS-rv-Vf . (10) 

By integrating this equation in velocity space, we obtain a steady-state continuity 
equation for the current density: 

V • J = 5 , (11) 

where J = / fv d 3 v and S = f S d 3 v. 

Multiplying Eq. (El by v and integrating again, we obtain 

J = J eq + tS - rV • n , (12) 

with II = n a/3 = f fv a vp d 3 v being the momentum-flux tensor and S = f Sv d 3 v 
the source current. The first term at the right-hand-side takes the form J eq = 
where u is the flow speed in the local equilibrium. For the case of a purely diffusive 
dynamics this term is zero. The second term contributes a shift tS to the flow speed, 
and it must also be set to zero for the case of pure diffusion. Finally, the pressure 
tensor shoud reduce to so that inserting (fT2l) into ELD. the diffusion equation 
is obtained. The first two conditions are automatically ensured by setting u = 0 in 
the local equilibrium. The third one, however, cannot be enforced exactly because 
of higher order contributions to the momemtum flux tensors. 

Indeed, by iterating the procedure, one obtains the general expression: 

OO 

^(-1)”t”V” +1 

n =0 

where n eq (”) = J f eq v ai ...v an d 3 v and S^ = J Sv ai ...v an d 3 1 >, are both tensors 
of order n. Note that for small Knudsen numbers (rV -C 1) higher-order terms 
become negligible and the series is convergent. We also observe that the solution is 
dictated by the equilibrium distribution and the source term, so that, by properly 
choosing both expressions, higher order terms can be canceled thus opening the way 
to high-accuracy solutions of the Poisson equation. 

In numerical practice, accuracy is typically improved by increasing the resolu¬ 
tion of the grid. This imposes correspondingly smaller time steps, which is clearly 
expensive, especially in three dimensions. High-order methods are meant to mit¬ 
igate the problem by achieving the same level of accuracy with many less grid 
points. They are typically based on higher-order stencils for the discretisation of 
the corresponding differential operators. 

In this work, we present a procedure whereby the same goal is achieved by 
drawing directly from an underlying lattice kinetic theory, i.e. by a suitable (non- 
Gaussian) extension of the local kinetic equilibria as combined with the use of larger 
sets of discrete speeds than those typically employed in standard Lattice Boltzmann 
theory. 


jjeq(n+l) r g(n+l) 


= S 


(13) 
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4. Kinetic approach to the Poisson equation 

In order for the Boltzmann kinetics to converge to the Poisson equation, the follow¬ 
ing conditions need to be met: 

n eq(l) = n eq(n>2) = q, ( 14 ) 

S (n >0 ) = q 


The equilibrium distribution can be expressed as a separable product of a function 
of the microscopic velocity (global uniform equilibrium) and a function of the time- 
spatial coordinates, namely: 

/ eq = (j)(y)<$>(x,t) , (15) 

and similarly for the source term: 

S = x(v)S(x,t) . (16) 


In the above, x(u) and 4>{v) are velocity dependent functions that need to be de¬ 
termined based on the conditions m above. To this purpose, we write: 


N 

(t>{v) =w{v)J2 b{ n )H n\v) , (17) 

n =0 

and, 

N 

X {v)=w{v)Y j ^ ) H^\v) , (18) 

n =0 


where b^ and c ^ are constant coefficients and TV is a cut-off truncating the Her- 
mite expansion. With reference to the one-dimensional case and an expansion up 
to seventh order (TV = 7), we obtain: 


(j>{v) ~ w{v) ( 1 — 


3c 2 — 6 c 2 v 2 + v 4 15c® — 45c 4 u 2 + 15c 2 u 4 — v 6 


, ( 19 ) 


and 


( 1 

v 2 

( 1+ 2 

1 - 2 

L c s J 


8c 4 24c® 

3c 2 — 6c 2 u 2 + v 4 15c® — A5c 4 v 2 + 15 c 2 v 4 — v 6 


8 c 4 


48 c? 


where the weight w(v ) for the one-dimensional case is defined as 


■j(v) = 


(2^c 2 )V2 


,~V 2 /1C 2 S 


( 20 ) 


( 21 ) 


With Eqs. (fl9l) and (12U1) . one can prove that the the following expressions for the 
moments are fulfilled: n eq (°) = 4>, II eq W = 0, n eq ( 2 ) = $c 2 , n eq ( rl>2 ) = 0, S® = S, 
g(n>o) _ q These expressions hold only up to TV = 7. 

To develop a lattice Boltzmann model and solve the Poisson equation numeri¬ 
cally, we need to impose a quadrature and find a corresponding set of finite velocity 
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vector s, Vi, such that the orthogonality conditions are fulfilled to the desired order 
| 24 | 25 | naine ly: 


w(v)H n (v)H m (v) dv = y WiH n (vi)H m (vi) 


N v 

2=0 


( 22 ) 


Here Wi are discrete weights and H n are the one-dimensional Hermite polynomials. 
The number of velocity vectors N v depends on the order of the expansion that we 
intend to reproduce with the quadrature. For expansions up to seventh order, we 
need at least 13 velocity vectors and weights. Each set of numbers also provides its 
own lattice sound speed c s , which goes from about 0.4 with N = 3 to nearly 1.2 
with N = 7. The details of these values are given in |Appendix A| 

Higher-order schemes are usually exposed to numerical instabilities, due to the 
appearance of additional modes in the dispersion relation. However, since our pro¬ 
cedure is designed to annihilate precisely these modes, we expect our approach to 
be stable. Furthermore, the method proposed here can also be used as a systematic 
technique of calculating weights and velocity vectors, such that one can achieve the 

desired level of accuracy for the computation of Laplacian operators, similar to the 
RT)T3T 1 

work proposed in Refs. | JU | ° ' . 

We have performed the theoretical analysis in one-dimension, however, exten¬ 
sions to two- and three-dimensions are straightforward, by using the tensorial form 
of the Hermite polynomials and by finding the discrete velocity vectors with the 
respective orthogonality conditions. For instance, in order to fulfill the conditions 
(HD, up to fifth order, the following expansions apply: 

15c! — 10 civ 2 + u 4 ' 


cf>(v) ~ w(v) ( 1 — 

and 

X(v) - w (v) 

where the weight w(v) is defined as 

w(v) = 


/ 1 

v 2 

1 + - 

3-^7 

V 2 

c 2 

L 4 J 


8 ci 


15c! - 10 dv 2 + v 4 


8cl 


(2^ C 2)3/2 


0 -v 2 /2c 


(23) 


(24) 


(25) 


For the three-dimensional case, we should calculate the corresponding set of discrete 
velocity vectors, u), such that the orthogonality condition is fulfilled up to fifth order, 
namely: 


N v 


w(v)H^\v)H^\v)d :i v = 

2=0 




(26) 


which is the three-dimensional version of Eq. (12211 . By solving these algebraic equa¬ 
tions, we find that we need at least 111 velocity vectors and 10 weights. Details are 
given in|Appcndix A[ 
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Fig. 1. Solution of the Poisson equation for different cut-off N and l = 1. The absolute error is 
computed A cp = (p n — (pt, where (p n is the potential calculated with the lattice Boltzmann model, 
and (pt is the analytical solution of Eq. @3- The system size is L = 128 lattice cells. 


5. Numerical Results 

As an application, we solve the Poisson equation 

V 2 $ = , (27) 

where p is the charge density and e is the electric permittivity. 

Replacing the equilibrium distribution / eq from Eq. (fl5l) into Eq. (fl3l) , we obtain: 

,S + c 2 tV 2 $~0 , (28) 

implying that S = pc 2 r/e. In the lattice Boltzmann model, one must take into 
account second order corrections to the relaxation time, r = 775 — 1 / 2 , where r u, the 
lattice Boltzmann relaxation time. 

From Eq. (fl3l) . it is appreciated that the spatial derivatives of the source also 
introduce errors in the solution. Thus, by annihilating higher order contributions, 
we are also eliminating spurious effects due to the derivatives of the source term. If 
the source term is smooth, this derivatives are negligible, however, for stiffer source 
terms, one can use the present approach to add accuracy to the solution. In order 
to investigate this issue, we solve the potential for the following charge density: 

p(x) = , (29) 

where the integer l controls the smoothness of the derivatives of the charge density, 
and L is the simulated length, with x £ [0, L). We use periodic boundary conditions 
for simplicity, and Tib = 1 (numerical units). 

We express the relative error as 

e = aL~ a (30) 

where the amplitude a and scaling exponent a (order of accuracy) both depend on 
N. From Fig.[TJ we see that for a relatively large system size, L = 128, the first order 
expansion is sufficient to produce satisfactory results, with errors around 0 . 02 % for 
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Fig. 2. Solution of the Poisson equation for different cut-off N and l = 4. The absolute error is 
computed A cj> = <f) n — <f>t, where ^> n is the potential calculated with the lattice Boltzmann model, 
and <j) t is the analytical solution of Eq. J27D . The system size is L = 128 lattice cells. 




Fig. 3. Solution of the Poisson equation for different cut-off N and l = 4. The absolute error is 
computed A</> = <j>n — <j>t, where </> n is the potential calculated with the lattice Boltzmann model, 
and 4>t is the analytical solution of Eq. m- The system size is L = 32 lattice cells. 


TV = 1 and TV = 3, and below 10 _3 % for TV = 5 and TV = 7. In Fig. [2l we see that 
by increasing the derivatives of the charge density (increasing l) 7 the error increases 
for all values of the cut-off TV. However, as we have mentioned before, the effect of 
removing higher-order errors becomes crucial for the case of small system sizes. By 
decreasing the number of lattice cells (see Fig. O, the discrepancies become visual 
and the errors are around 5% for TV = 1,3, 1% for TV = 5, and less than 0.4% for 
TV = 7. 

Finally, by further reducing the lattice resolution, L = 8 cells, and considering 
only one oscillation, 1 = 1. we see that the errors for the first order expansion, 
TV = 1, are around 7%, while for the highest values of TV = 7 they remain below 
0.5%. The fact that one can obtain very small errors by increasing the order of 
the expansion in the velocity space of the equilibrium distribution and source term, 
opens up the possibility of highly accurate solutions on very small grids. 

Increasing the order of the expansion inevitably implies an overhead of numerical 
operations per time step, thereby slowing down the simulations. On the other hand, 
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Fig. 4. Solution of the Poisson equation for different cut-off N and 1 = 1. The absolute error is 
computed A </> = <f) n — <f>t, where ^> n is the potential calculated with the lattice Boltzmann model, 
and 4>t is the analytical solution of Eq. @3- The system size is L = 8 lattice cells. 



L (cells) 



Fig. 5. (left) Relative error as a function of the system size L for each cut-off value N of the 
Hermite polynomials expansion. The lines denote the fitting curves using the general expression, 
e = clnL~ a (right). Number of iterations as a function of the system size L for each cut-off value 
of N. The lines denote the fitting curve using the general expression, M = bj^L 2 , M being the 
number of time-steps. 


higher order lattices also bring aditional side-benefits such as a higher sound speed 
c s , leading to a larger diffusion coefficient and, consequently, to a faster path to the 
steady-state solution. In order to highlight the concrete advantages of our approach, 
we next implement the same simulations for different orders N and inspect the CPU 
time, number of iterations, and system size, required to achieve a specific level of 
accuracy. 


5.1. Computational performance 

Let us next change the order of the expansion and analyse how the relative error 
decreases with the system size. From Fig. [5] (left), we see that, according to our 
expectations, by increasing the order of the expansion, the order of convergence, 
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Table 1. Parameters ajy and bjy for each value of N. Also 
shown are the accuracy exponent a .jy and the number of 
lattice sites updated per second, £jv, in Msites/s (Million 
sites per second). 


N 

a N 

b N 

£jv (Msites/s) 

Ijv/^jv 

&N 

1 

1.645 

10.98 

26.2 

2.39 

1 

3 

1.645 

10.98 

21.4 

1.95 

1 

5 

33.11 

5.46 

5.44 

0.997 

3 

7 

992.3 

3.23 

5.76 

1.78 

5 


measured by the exponent a, also increases and in a very substantial way. From 
the same figure (right), we can also observe that for all values of TV , the number 
of time-steps M scales quadratically with the lattice size, and decreases with the 
order of the expansion TV. The former feature is the typical signature of diffusion 
process, while the second reflects the fact that higher order lattice deliver a larger 
sound speed, hence they take less time-steps to achieve a given level of accuracy. 

The computational time required to achieve a given accuracy clearly depends on 
the computational speed of the model on each given machine. This is conveniently 
expressed in terms of the number of lattice sites updates per CPU second, £jv- This 
quantity is expected to be a decreasing function of TV, since higher-order quadratures 
require more operations per lattice site. 

Actual measurements give the values reported in Table [T] Although no clear 
trend can be extracted from these measurements, it is observed that the cases 
TV = 5 and TV = 7 are about five times more expensive than the cases TV = 1 and 
TV = 3. Also to be noted that £7 > £5, which is counter-intuitive. However, upon 
inspecting the velocity vectors (see Appendix | Appendix A| ), one can appreciate that 
for TV = 7 the velocities are consecutive (no velocity gap), while for TV = 5 they 
are not. Since, memory access is faster for consecutive elements, we expect the case 
TV = 7 to perform better than TV = 5. More generally, the dependence on TV appears 
to be highly dependent on memory access patterns, which are not easily quantified 
through simple scaling relations. 

The total CPU time, Tl.at, it takes to attain steady-state for a given system size 
L can be written as follows: 

tL,N = 7 —L 3 . (31) 

Uv 

In the above, the cubic dependence derives from the L 2 diffusive scaling, as combined 
with the number of sites L (the space-time computational volume of a diffusive 
process scales like L 2+D in D spatial dimensions). Finally, is a relative scale for 
the number of time-steps, which is expected to decrease at increasing TV because 
of the increasing sound speed. The actual values of obtained in the simulations 
are reported in Table [T] As an example, for the case L = 100, we find that the 
simulations take 0.42 (65629 time steps), 0.51 (65629 time steps), 1.0 (33769 time 
steps), and 0.56 (20079 time steps) seconds, for TV = 1,3, 5,7 respectively. 

This shows that the total computational time remains within a factor two across 
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Fig. 6. Computational time tjy and system size L needed to recover the solution of the Poisson 
equation with a relative error of 0.1%, for different N. 


all values of TV , while the error with respect to the analytical solution decreases 
dramatically at increasing TV (see Fig. 0. Hence, we can conclude that increasing 
TV leads to a dramatic boost of accuracy at a very moderate extra computational 
cost. Finally, we have inspected the computational time to steady-state for a given 
accuracy, at changing the size of the problem and the order of the quadrature. By 
combining the previous relations, m and (ED, we obtain: 

M«) = ■ < 32 > 

The plot reports the computational time and system size needed to achieve a pre¬ 
scribed relative error, for the case in point e = 10 -3 . From fig. [Gj we see that by 
increasing the order of the expansion in Hermite polynomials, it is possible to reduce 
the computational time by several orders of magnitude (up to six). 

These data refer specifically to the solution of the Poisson equation, but we have 
reasons to believe that similar conclusions would hold for other partial differential 
equations with a kinetic-theory background, including many linear non-linear partial 
differential equations of great relevance in many branches of modern science. 

5.1.1. Comparison with multi-grid, methods 

In order to gain information about the performance of this approach, as compared 
with existing tools to solve the Pois son eq uation, we have implemented simulations 
using the multi-grid (MG) method |32 | 33| _ num b er of nodes in the MG simu¬ 
lations is denoted by L, like for the lattice Boltzmann method. Fixing l = 1 (from 
Eq. (l29l) b in Table [2] we show the relative errors per node, e/L, together with the 
corresponding computational time. Note that for IV = 3, the MG method is faster 
than our approach, with the same second order accuracy. By increasing TV, however, 
the MG becomes increasingly faster, but also increasingly less accurate at a given 
grid size L. For instance, in order to attain the same accuracy of the case TV = 7 on 
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Table 2. Average relative error per node, e jL. for the solution of the Pois¬ 
son equation using different orders N and the multi-grid method {MG). 
Inside the parenthesis, we report the computational time in milliseconds. 


L 

N = 3 

N = 5 

N = 7 

MG 

16 

10“ 2 (2) 

9 x 10“ 4 (2) 

10" 4 (2) 

icr 2 (< i) 

32 

3 x 10 -a (4) 

6 x 10 —5 (4) 

2 x 10 -h (6 ) 

3 x 10 3 (< 1) 

64 

8 x 10 4 (26) 

4 x 10 t> (29) 

3 x 10 ts (38) 

8 x icr 4 (< i) 

2 14 

- 

- 

- 

2 x 10“ s (12) 


L = 64, the MG solver requires L = 16384 grid points, being still faster by about 
a factor 3. Given that N = 7 requires 13 arrays, the LB memory saving is about 
a factor 20. These figures indicate that the present high-order LB Poisson solver is 
consistently slower than MG, however it also takes correspondingly less memory to 
achieve high levels of accuracy. Given that MG represents the state-of-the art for 
fast Poisson solvers, these can be regarded as satisfactory results. 

In summary, we have shown that the present LB approach provides a viable 
option for the solution of Poisson equation, in the sense of simplicity and possibly 
also in terms of memory demand at high accuracy. 


6. Conclusions 

Summarizing, we have studied the effect of higher order terms on the steady-state 
solution of the diffusion equation, using a lattice version of the Boltzmann equa¬ 
tion with high-order (non-gaussian) equilibria. For the numerical solution, we have 
developed a higher order lattice Boltzmann model capable of reproducing up to 
seventh order moment of the equilibrium distribution and source charge. This per¬ 
mits to cancel exactly error terms up to seventh order. For the validation and study 
of our approach, we have solved the one-dimensional Poisson equation, and found 
that the computational time to steady-state, at a prescribed level of accuracy, can 
be reduced by up to six orders of magnitude as compared to standard LB formu¬ 
lations. The six orders gain in computational time is due to the fact that higher 
order lattices permit to utilize much smaller grid sizes. Moreover, since they contain 
a large number of discrete velocities, they allow faster propagation within the lat¬ 
tice, hence a correspondingly faster attainement of the steady-state. We have also 
compared with state-of-the art multigrid solvers, and found that high-order LB is 
nearly competitive in efficiency at a lower memory demand. 

Extensions to two and three-dimensions, as well as time-dependent diffusion and 
advection-diffusion equations will make the object of future research. 
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Appendix A. Quadrature of the Boltzmann Equation 

Using Eq. (l22l) for each order of expansion A, we can calculate the velocity vectors 
Vi with the respective weights Wj. In addition, each set of Wi and provides a 
characteristic normalised speed c s , c 2 = J2iLo w i v i- Here we only show the cases 
for N = 3, 5, and 7, since the simulations for N = 1 were performed with the same 
lattice than for N = 3. For N = 3, we have 


w 0 = 0.6366469031260781628443461 

Wl = 0.18141458774368577505004149 

w 2 = 0.18141458774368577505004149 

w 3 = 0.0002619606932751435277854615 (A.l) 

w 4 = 0.0002619606932751435277854615 

v 0 = 0, vi = 1, v 2 = -1, v 3 = 3, v 4 = -3 

c 2 s = l-y/2jl , 


for N = 5, 


w 0 = 0.45813515550767658573 

wi = 0.23734280857794043891 

w 2 = 0.23734280857794043891 

w 3 = 0.032324653788654934092 

w 4 = 0.032324653788654934092 

w 5 = 0.0012640621515365148385 

w 6 = 0.0012640621515365148385 

w 7 = 8.977280298192933351 x 10“ 7 

w 8 = 8.977280298192933351 x 10 -7 

v 0 = 0, wl = 1, v 2 = -1, v 3 = 2, v 4 = -2, 

v 5 = 3, v 6 = -3, v 7 = 5, v 8 = -5 

c 2 s = 0.86952909818721338884964163917209 


(A.2) 


I 


I 
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Table 3. Velocity vectors Vi with their corresponding weights Wi 
for N = 5. There are 111 velocity vectors and 10 different weights. 


Wi 

Vi 

0.15014405 

(0,0,0) 

0.02500399 

(±1,0, 0),(0, ±1,0), (0, 0,±1) 

0.04505812 

(±1, ±1,0), (0, ±1, ±1), (±1, 0, ±1) 

7.81706951 x 10 7 

(±3, ±3, ±3) 

9.79290909 x 10“® 

(±3, ±3, 0), (±3, 0, ±3), (0, ±3, ±3) 

1.81207346 x 10 -4 

(±3, ±1,0), (±3, 0, ±1), (0, ±1, ±3) 

(±1, ±3,0), (±1, 0, ±3), (0, ±3, ±1) 

6.03109965 x 10“ 4 

(±2, ±2, 0), (±2, 0, ±2), (0, ±2, ±2) 

3.49417975 x 10 -a 

(±2, ±1, ±1), (±1, ±1, ±2), (±1, ±2, ±1) 

1.05490305 x 10'^ 

(±2, 0, 0), (0, ±2,0), (0, 0, ±2) 

4.50016852 x 10“ 5 

(±3,0,0), (0, ±3,0), (0,0, ±3) 


and finally, for N = 7, 

w 0 = 0.3455934552621565 
Wl = 0.2374599218260301 
w 2 = 0.2374599218260301 
w 3 = 0.07705730993964580 
w 4 = 0.07705730993964580 
w 5 = 0.011801423732312036 
w 6 = 0.011801423732312036 
w 7 = 0.0008552466009513439 

w 8 = 0.0008552466009513439 (A.3) 

w 9 = 0.00002884934614927074 
w w = 0.00002884934614927074 
wh = 5.209238332209471 x 1CT 7 
w 12 = 5.209238332209471 x 1(T 7 
v 0 =0,vi = 1, v 2 = -1, V 3 = 2, v 4 = -2 
v 5 = 3, ^6 = -3, wj = 4, v 8 = -4, v 9 = 5 
v\o = -5, vu = 6, V 12 = -6 
c 2 s = 1.1544053947399681272395977588380 . 

For the three dimensional case, up to fifth order, we can find the velocity vectors 
Vi and weights Wi by solving Eq. (1261) . The values are shown in Tabic [3] 
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